library(tidyverse)
library(sensemakr)
library(estimatr)
library(texreg)
library(corrr)
library(lmerTest)
library(lme4)
library(interplot)
library(margins)
library(stargazer)
library(patchwork)
library(countrycode)
library(ggrepel)
library(regweight)

## Functions

zero1 <- function(x, minx=NA, maxx=NA){
  res <- NA
  if(is.na(minx)) res <- (x - min(x,na.rm=T))/(max(x,na.rm=T) -min(x,na.rm=T))
  if(!is.na(minx)) res <- (x - minx)/(maxx -minx)
  res
}



load("cleaned_data/study2_data.RData")

cultmodels <- df_na %>% drop_na(cultatts,kinship_score,polengage,X001, X002 ,ln_time_obs_ea ,small_scale,agriculture ,malariaindex,S002,S003)

culture1 <- lm_robust(cultatts ~ kinship_score,
                      fixed_effects = ~ as.factor(S002) + as.factor(S003),
                      data = cultmodels,
                      clusters = group,
                      weights = S017,
                      se_type = "stata")
culture2 <- lm_robust(cultatts ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                      fixed_effects = ~ as.factor(S002) + as.factor(S003),
                      data = cultmodels,
                      clusters = group,
                      weights = S017,
                      se_type = "stata")

culture3 <- lm_robust(cultatts ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex + as.factor(F025) + X025R,
                      fixed_effects = ~ as.factor(S002) + as.factor(S003),
                      data = cultmodels,
                      clusters = group,
                      weights = S017,
                      se_type = "stata")


culture4 <- lm_robust(cultatts ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex + as.factor(F025) + X025R + authoritarianism1,
                      fixed_effects = ~ as.factor(S002) + as.factor(S003),
                      data = cultmodels,
                      clusters = group,
                      weights = S017,
                      se_type = "stata")

culture5 <- lm_robust(cultatts ~ kinship_score * polengage + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                      fixed_effects = ~ as.factor(S002) + as.factor(S003),
                      data = cultmodels,
                      clusters = group,
                      weights = S017,
                      se_type = "stata")

texreg(list(culture1, culture2, culture3, culture5), omit.coef=c('X001|X002|ln_time_obs_ea|small_scale|malariaindex|F025|agriculture|X025R'),
       custom.gof.rows = list("Country fixed effects" = c("Yes","Yes","Yes","Yes"),
                              "Wave fixed effects" = c("Yes","Yes","Yes","Yes"),
                              "Individual-level controls" = c("No","Yes","Yes","Yes"),
                              "Ethnicity-level controls" = c("No","Yes","Yes","Yes"),
                              "Religion-level controls" = c("No","No", "Yes","No")),
       custom.coef.names = c("Kinship Score", "Political Engagement","Kinship * Pol Engagement"), include.adjrs = FALSE,
       include.rmse = FALSE,
       scalebox = .9,include.ci=F,file = "output/wvs_cultural_final_table.tex",
       override.pvalues = list(culture1$p.value/2,culture2$p.value/2,culture3$p.value/2,culture5$p.value/2),
       caption = "Cultural Models",use.packages = F)

econattsmodels <- df_na %>% drop_na(econatts,kinship_score,polengage,X001, X002 ,ln_time_obs_ea ,small_scale,agriculture ,malariaindex,S002,S003)


econ1 <- lm_robust(econatts ~ kinship_score,
                   fixed_effects = ~ as.factor(S002) + as.factor(S003),
                   data = econattsmodels,
                   clusters = group,
                   weights = S017,
                   se_type = "stata")

econ2 <- lm_robust(econatts ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                   fixed_effects = ~ as.factor(S002) + as.factor(S003),
                   data = econattsmodels,
                   clusters = group,
                   weights = S017,
                   se_type = "stata")

econ3 <- lm_robust(econatts ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex + as.factor(F025) + X025R,
                   fixed_effects = ~ as.factor(S002) + as.factor(S003),
                   data = econattsmodels,
                   clusters = group,
                   weights = S017,
                   se_type = "stata")

econ4 <- lm_robust(econatts ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex + as.factor(F025) + authoritarianism1,
                   fixed_effects = ~ as.factor(S002) + as.factor(S003),
                   data = econattsmodels,
                   clusters = group,
                   weights = S017,
                   se_type = "stata")

econ5 <- lm_robust(econatts ~ kinship_score * polengage + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                   fixed_effects = ~ as.factor(S002) + as.factor(S003),
                   data = econattsmodels,
                   clusters = group,
                   weights = S017,
                   se_type = "stata")

texreg(list(econ1, econ2, econ3, econ5), omit.coef=c('X001|X002|ln_time_obs_ea|small_scale|malariaindex|F025|agriculture|X025R'),
       custom.gof.rows = list("Country fixed effects" = c("Yes","Yes","Yes","Yes"),
                              "Wave fixed effects" = c("Yes","Yes","Yes","Yes"),
                              "Individual-level controls" = c("No","Yes","Yes","Yes"),
                              "Ethnicity-level controls" = c("No","Yes","Yes","Yes"),
                              "Religion-level controls" = c("No","No", "Yes","No")),
       custom.coef.names = c("Kinship Score", "Political Engagement","Kinship * Pol Engagement"), include.adjrs = FALSE,
       include.rmse = FALSE,
       scalebox = .9,include.ci=F,file = "output/wvs_econ_final_table.tex",
       override.pvalues = list(econ1$p.value/2,econ2$p.value/2,econ3$p.value/2,econ5$p.value/2),
       caption = "Economic Models",use.packages = F)

## sensitivity analysis

df_ses <- df_na %>% mutate(Muslim = ifelse(F025 == "Muslim", 1, 0))

culture3_sens <- lm(cultatts ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex + Muslim,
                    data = df_ses,
                    weights = S017)

sensitivity <- sensemakr::sensemakr(model = culture3_sens, 
                                    treatment = "kinship_score",
                                    benchmark_covariates = "Muslim",
                                    kd = 1:4,
                                    q = 1,
                                    alpha = 0.05, 
                                    reduce = TRUE)
summary(sensitivity)
plot(sensitivity)
cor(df_ses$kinship_score, df_ses$Muslim)



# Number of ethnicities in each country -----------------------------------

group_unique <- df_na %>% dplyr::select(group, S003) %>% 
  group_by(S003,group) %>% summarise(n = n()) %>%   mutate(Freq = n/sum(n))

group2 <- df_na %>% dplyr::select(S003) %>% 
  group_by(S003) %>% summarise(n = n()) %>%   mutate(Freq = n/sum(n))

unique(group_unique$S003)

group_unique1 <- table(group_unique$S003)
as.data.frame(group_unique1) %>% filter(Freq!=1) -> tab
#df_na %>% filter(S003 %in% as.numeric(as.character(tab$Var1))) -> out
#hist(out$kinship_score)
#out %>% group_by(S003)


# Models for each attitude question ---------------------------------------

homosexuality1 <- lm_robust(homosexualitynotok ~ kinship_score,
                            fixed_effects = ~ as.factor(S002) + as.factor(S003),
                            data = df_na,
                            clusters = group,
                            weights = S017,
                            se_type = "stata")

homosexuality2 <- lm_robust(homosexualitynotok ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                            fixed_effects = ~ as.factor(S002) + as.factor(S003),
                            data = df_na,
                            clusters = group,
                            weights = S017,
                            se_type = "stata")

homosexuality3 <- lm_robust(homosexualitynotok ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex + as.factor(F025) + X025R,
                            fixed_effects = ~ as.factor(S002) + as.factor(S003),
                            data = df_na,
                            clusters = group,
                            weights = S017,
                            se_type = "stata")

homosexuality5 <- lm_robust(homosexualitynotok ~ kinship_score * polengage + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                            fixed_effects = ~ as.factor(S002) + as.factor(S003),
                            data = df_na,
                            clusters = group,
                            weights = S017,
                            se_type = "stata")

texreg(list(homosexuality1, homosexuality2, homosexuality3, homosexuality5), omit.coef=c('X001|X002|ln_time_obs_ea|small_scale|malariaindex|F025|agriculture|X025R'),
       custom.gof.rows = list("Country fixed effects" = c("Yes","Yes","Yes","Yes"),
                              "Wave fixed effects" = c("Yes","Yes","Yes","Yes"),
                              "Individual-level controls" = c("No","Yes","Yes","Yes"),
                              "Ethnicity-level controls" = c("No","Yes","Yes","Yes"),
                              "Religion-level controls" = c("No","No", "Yes","No")),
       custom.coef.names = c("Kinship Score", "Political Engagement","Kinship * Pol Engagement"), include.adjrs = FALSE,
       include.rmse = FALSE,
       scalebox = .9,include.ci=F,file = "output/wvs_homo_final_table.tex",
       override.pvalues = list(homosexuality1$p.value/2,homosexuality2$p.value/2,homosexuality3$p.value/2,homosexuality5$p.value/2),
       caption = "Gay Attitude Models",use.packages = F)




abort1 <- lm_robust(abortionnotok ~ kinship_score,
                    fixed_effects = ~ as.factor(S002) + as.factor(S003),
                    data = df_na,
                    clusters = group,
                    weights = S017,
                    se_type = "stata")

abort2 <- lm_robust(abortionnotok ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                    fixed_effects = ~ as.factor(S002) + as.factor(S003),
                    data = df_na,
                    clusters = group,
                    weights = S017,
                    se_type = "stata")

abort3 <- lm_robust(abortionnotok ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex + as.factor(F025) + X025R,
                    fixed_effects = ~ as.factor(S002) + as.factor(S003),
                    data = df_na,
                    clusters = group,
                    weights = S017,
                    se_type = "stata")

abort5 <- lm_robust(abortionnotok ~ kinship_score * polengage + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                    fixed_effects = ~ as.factor(S002) + as.factor(S003),
                    data = df_na,
                    clusters = group,
                    weights = S017,
                    se_type = "stata")

texreg(list(abort1, abort2, abort3, abort5), omit.coef=c('X001|X002|ln_time_obs_ea|small_scale|malariaindex|F025|agriculture|X025R'),
       custom.gof.rows = list("Country fixed effects" = c("Yes","Yes","Yes","Yes"),
                              "Wave fixed effects" = c("Yes","Yes","Yes","Yes"),
                              "Individual-level controls" = c("No","Yes","Yes","Yes"),
                              "Ethnicity-level controls" = c("No","Yes","Yes","Yes"),
                              "Religion-level controls" = c("No","No", "Yes","No")),
       custom.coef.names = c("Kinship Score", "Political Engagement","Kinship * Pol Engagement"), include.adjrs = FALSE,
       include.rmse = FALSE,
       scalebox = .9,include.ci=F,file = "output/wvs_abort_final_table.tex",
       override.pvalues = list(abort1$p.value/2,abort2$p.value/2,abort3$p.value/2,abort5$p.value/2),
       caption = "Abortion Models",use.packages = F)


respon1 <- lm_robust(econresponsibility ~ kinship_score,
                     fixed_effects = ~ as.factor(S002) + as.factor(S003),
                     data = df_na,
                     clusters = group,
                     weights = S017,
                     se_type = "stata")

respon2 <- lm_robust(econresponsibility ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                     fixed_effects = ~ as.factor(S002) + as.factor(S003),
                     data = df_na,
                     clusters = group,
                     weights = S017,
                     se_type = "stata")

respon3 <- lm_robust(econresponsibility ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex + as.factor(F025) + X025R,
                     fixed_effects = ~ as.factor(S002) + as.factor(S003),
                     data = df_na,
                     clusters = group,
                     weights = S017,
                     se_type = "stata")

respon5 <- lm_robust(econresponsibility ~ kinship_score * polengage + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                     fixed_effects = ~ as.factor(S002) + as.factor(S003),
                     data = df_na,
                     clusters = group,
                     weights = S017,
                     se_type = "stata")

texreg(list(respon1, respon2, respon3, respon5), omit.coef=c('X001|X002|ln_time_obs_ea|small_scale|malariaindex|F025|agriculture|X025R'),
       custom.gof.rows = list("Country fixed effects" = c("Yes","Yes","Yes","Yes"),
                              "Wave fixed effects" = c("Yes","Yes","Yes","Yes"),
                              "Individual-level controls" = c("No","Yes","Yes","Yes"),
                              "Ethnicity-level controls" = c("No","Yes","Yes","Yes"),
                              "Religion-level controls" = c("No","No", "Yes","No")),
       custom.coef.names = c("Kinship Score", "Political Engagement","Kinship * Pol Engagement"), include.adjrs = FALSE,
       include.rmse = FALSE,
       scalebox = .9,include.ci=F,file = "output/wvs_responsibility_final_table.tex",
       override.pvalues = list(respon1$p.value/2,respon2$p.value/2,respon3$p.value/2,respon5$p.value/2),
       caption = "Econ Responsibility Models",use.packages = F)


income1 <- lm_robust(incomediffs ~ kinship_score,
                     fixed_effects = ~ as.factor(S002) + as.factor(S003),
                     data = df_na,
                     clusters = group,
                     weights = S017,
                     se_type = "stata")

income2 <- lm_robust(incomediffs ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                     fixed_effects = ~ as.factor(S002) + as.factor(S003),
                     data = df_na,
                     clusters = group,
                     weights = S017,
                     se_type = "stata")

income3 <- lm_robust(incomediffs ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex + as.factor(F025) + X025R,
                     fixed_effects = ~ as.factor(S002) + as.factor(S003),
                     data = df_na,
                     clusters = group,
                     weights = S017,
                     se_type = "stata")

income5 <- lm_robust(incomediffs ~ kinship_score * polengage + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                     fixed_effects = ~ as.factor(S002) + as.factor(S003),
                     data = df_na,
                     clusters = group,
                     weights = S017,
                     se_type = "stata")

texreg(list(income1, income2, income3, income5), omit.coef=c('X001|X002|ln_time_obs_ea|small_scale|malariaindex|F025|agriculture|X025R'),
       custom.gof.rows = list("Country fixed effects" = c("Yes","Yes","Yes","Yes"),
                              "Wave fixed effects" = c("Yes","Yes","Yes","Yes"),
                              "Individual-level controls" = c("No","Yes","Yes","Yes"),
                              "Ethnicity-level controls" = c("No","Yes","Yes","Yes"),
                              "Religion-level controls" = c("No","No", "Yes","No")),
       custom.coef.names = c("Kinship Score", "Political Engagement","Kinship * Pol Engagement"), include.adjrs = FALSE,
       include.rmse = FALSE,
       scalebox = .9,include.ci=F,file = "output/wvs_income_final_table.tex",
       override.pvalues = list(income1$p.value/2,income2$p.value/2,income3$p.value/2,income5$p.value/2),
       caption = "Income Difference Models",use.packages = F)



# Histogram ---------------------------------------------------------------

ggplot(df_na, aes(x=kinship_score)) + 
  geom_histogram(binwidth = .1)+
  theme_bw()+
  ylab("Count")+
  xlab("Kinship Tightness")

ggsave("output/kinship_wvs.png", width = 5, height = 4, dpi = 600, device = 'png')



# Coefficient Plot --------------------------------------------------------

c1 <- broom::tidy(culture1, conf.int = TRUE)
c1$dv <- "Culture"
c1$model <- "Kinship Tightness"
c1$order <- "1"

c2 <- broom::tidy(culture2, conf.int = TRUE)
c2$dv <- "Culture"
c2$model <- "Kinship Tightness"
c2 <- c2[1,]
c2$order <- "2"

c3 <- broom::tidy(culture3, conf.int = TRUE)
c3$dv <- "Culture"
c3$model <- "Extended Models"
c3 <- c3[1,]
c3$order <- "3"

c4 <- broom::tidy(culture5, conf.int = TRUE)
c4$dv <- "Culture"
c4$model <- "Extended Models"
c4 <- c4[1,]
c4$order <- "4"



e1 <- broom::tidy(econ1, conf.int = TRUE)
e1$dv <- "Economic"
e1$model <- "Kinship Tightness"
e1$order <- "1"

e2 <- broom::tidy(econ2, conf.int = TRUE)
e2$dv <- "Economic"
e2$model <- "Kinship Tightness"
e2 <- e2[1,]
e2$order <- "2"

e3 <- broom::tidy(econ3, conf.int = TRUE)
e3$dv <- "Economic"
e3$model <- "Extended Models"
e3 <- e3[1,]
e3$order <- "3"

e4 <- broom::tidy(econ4, conf.int = TRUE)
e4$dv <- "Economic"
e4$model <- "Extended Models"
e4 <- e4[1,]
e4$order <- "4"


total <- rbind(c1, c2, c3, e1, e2, e3)

total <- total %>% mutate(order = order,
                          model = fct_rev(model))

ggplot(total, aes(order, estimate)) +
  geom_point(aes(), position = position_dodge(width = .55), size = 2)+
  geom_pointrange(aes(ymin = conf.low, ymax = conf.high), alpha = 0.9,
                  position = position_dodge(width = .55))+
  theme_bw()+
  facet_wrap(~dv, nrow = 1)+
  geom_hline(yintercept = 0, colour = "grey60", linetype = 2)+
  ylab("Estimate")+
  xlab("Specification")

ggsave("output/coefficient_plot3.png", units = "in", width = 7, height = 3.5, dpi = 600, device = 'png')




# Schwartz Values ---------------------------------------------------------


cultmodels_schwartz <- df_na %>% drop_na(cultatts,econatts,kinship_score,polengage,X001, 
                                         X002 ,ln_time_obs_ea ,small_scale,agriculture ,
                                         malariaindex,S002,S003,opennesstochange,conservation,selfenhancement,selftranscendence)


m1a <- lm_robust(opennesstochange ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                 fixed_effects = ~ as.factor(S002) + as.factor(S003),
                 data = cultmodels_schwartz,
                 clusters = group,
                 weights = S017,
                 se_type = "stata")
summary(m1a)

m1b <- lm_robust(conservation ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                 fixed_effects = ~ as.factor(S002) + as.factor(S003),
                 data = cultmodels_schwartz,
                 clusters = group,
                 weights = S017,
                 se_type = "stata")
summary(m1b)

m1c <- lm_robust(selfenhancement ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                 fixed_effects = ~ as.factor(S002) + as.factor(S003),
                 data = cultmodels_schwartz,
                 clusters = group,
                 weights = S017,
                 se_type = "stata")
summary(m1c)

m1d <- lm_robust(selftranscendence ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                 fixed_effects = ~ as.factor(S002) + as.factor(S003),
                 data = cultmodels_schwartz,
                 clusters = group,
                 weights = S017,
                 se_type = "stata")
summary(m1d)


texreg(list(m1a, m1b, m1c, m1d), custom.coef.names = c("Kinship Score","Gender", "Age", "Date entered in Ethnographic Atlas","Small Scale Group","Agriculture", "Malaria Index"), custom.model.names=c("Openness to Change", "Conservation", "Self-Enhancement", "Self-Transcendence"), file = "output/wvs_moderation.tex",custom.note = "All models include country-of-residence and wave fixed effects. SEs are clustered by ethnic group.",caption = "Kinship Predicting Moral Values",scalebox = .9,include.ci=F,override.pvalues = list(m1a$p.value/2,m1b$p.value/2,m1c$p.value/2,m1d$p.value/2),table = F,use.packages = F)


m2a  <- lm_robust(cultatts ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                  fixed_effects = ~ as.factor(S002) + as.factor(S003),
                  data = cultmodels_schwartz,
                  clusters = group,
                  weights = S017,
                  se_type = "stata")
summary(m2a)

m2b  <- lm_robust(cultatts ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex +
                    opennesstochange + conservation + selfenhancement + selftranscendence,
                  fixed_effects = ~ as.factor(S002) + as.factor(S003),
                  data = cultmodels_schwartz,
                  clusters = group,
                  weights = S017,
                  se_type = "stata")
summary(m2b)


m2c  <- lm_robust(econatts ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                  fixed_effects = ~ as.factor(S002) + as.factor(S003),
                  data = cultmodels_schwartz,
                  clusters = group,
                  weights = S017,
                  se_type = "stata")

m2d  <- lm_robust(econatts ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex +
                    opennesstochange + conservation + selfenhancement + selftranscendence,
                  fixed_effects = ~ as.factor(S002) + as.factor(S003),
                  data = cultmodels_schwartz,
                  clusters = group,
                  weights = S017,
                  se_type = "stata")

summary(m2b)


texreg::texreg(list(m2a, m2b, m2c, m2d), custom.coef.names = c("Kinship Score","Gender", "Age","Date entered in Ethnographic Atlas","Small Scale Group","Agriculture", "Malaria Index", "Openness to Change", "Conservation", "Self-Enhancement", "Self-Transcendence"), file = "output/wvs_moderation2.tex",custom.note = "All models include country-of-residence and wave fixed effects. SEs are clustered by ethnic group.",caption = "Kinship Models with Moral Values",scalebox = .9,include.ci=F,override.pvalues = list(m2a$p.value/2,m2b$p.value/2,m2c$p.value/2,m2d$p.value/2),table = F,use.packages = F)



cultmodels_schwartz <- df_na %>% drop_na(cultatts,econatts,kinship_score,polengage,X001, X002 ,ln_time_obs_ea ,small_scale,agriculture ,malariaindex,S002,S003,authoritarianism1)


n1a <- lm_robust(authoritarianism1 ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                 fixed_effects = ~ as.factor(S002) + as.factor(S003),
                 data = cultmodels_schwartz,
                 clusters = group,
                 weights = S017,
                 se_type = "stata")

summary(n1a)

n1b <- lm_robust(cultatts ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                 fixed_effects = ~ as.factor(S002) + as.factor(S003),
                 data = cultmodels_schwartz,
                 clusters = group,
                 weights = S017,
                 se_type = "stata")

n1c <- lm_robust(cultatts ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex + authoritarianism1,
                 fixed_effects = ~ as.factor(S002) + as.factor(S003),
                 data = cultmodels_schwartz,
                 clusters = group,
                 weights = S017,
                 se_type = "stata")

n1d <- lm_robust(econatts ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                 fixed_effects = ~ as.factor(S002) + as.factor(S003),
                 data = cultmodels_schwartz,
                 clusters = group,
                 weights = S017,
                 se_type = "stata")

n1e <- lm_robust(econatts ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex + authoritarianism1,
                 fixed_effects = ~ as.factor(S002) + as.factor(S003),
                 data = cultmodels_schwartz,
                 clusters = group,
                 weights = S017,
                 se_type = "stata")

summary(n1b)
df_na$trust
texreg::texreg(list(n1a, n1b,n1c,n1d,n1e), custom.coef.names = c("Kinship Score","Gender", "Age","Date entered in Ethnographic Atlas","Small Scale Group","Agriculture", "Malaria Index", "Authoritarianism"), file = "output/wvs_authoritarianism.tex",custom.note = "All models include country-of-residence and wave fixed effects. SEs are clustered by ethnic group.",caption = "Kinship Models with Authoritarianism",scalebox = .9,include.ci=F,override.pvalues = list(n1a$p.value/2,n1b$p.value/2, n1c$p.value/2, n1d$p.value/2, n1e$p.value/2),
               table = F,use.packages = F,
               custom.model.names = c("Authoritarianism", "Cultural Attitudes",  "Cultural Attitudes",  "Economic Attitudes",  "Economic Attitudes"))


cultmodels_schwartz <- df_na %>% drop_na(cultatts,econatts,kinship_score,polengage,X001, X002 ,ln_time_obs_ea ,small_scale,agriculture ,malariaindex,S002,S003,trust)


n1a <- lm_robust(trust ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                 fixed_effects = ~ as.factor(S002) + as.factor(S003),
                 data = cultmodels_schwartz,
                 clusters = group,
                 weights = S017,
                 se_type = "stata")

summary(n1a)

n1b <- lm_robust(cultatts ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                 fixed_effects = ~ as.factor(S002) + as.factor(S003),
                 data = cultmodels_schwartz,
                 clusters = group,
                 weights = S017,
                 se_type = "stata")

n1c <- lm_robust(cultatts ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex + trust,
                 fixed_effects = ~ as.factor(S002) + as.factor(S003),
                 data = cultmodels_schwartz,
                 clusters = group,
                 weights = S017,
                 se_type = "stata")

n1d <- lm_robust(econatts ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex,
                 fixed_effects = ~ as.factor(S002) + as.factor(S003),
                 data = cultmodels_schwartz,
                 clusters = group,
                 weights = S017,
                 se_type = "stata")

n1e <- lm_robust(econatts ~ kinship_score + X001 + X002 + ln_time_obs_ea + small_scale + agriculture + malariaindex + trust,
                 fixed_effects = ~ as.factor(S002) + as.factor(S003),
                 data = cultmodels_schwartz,
                 clusters = group,
                 weights = S017,
                 se_type = "stata")

texreg::texreg(list(n1a, n1b,n1c,n1d,n1e), custom.coef.names = c("Kinship Score","Gender", "Age","Date entered in Ethnographic Atlas","Small Scale Group","Agriculture", "Malaria Index", "Trust"), file = "output/wvs_trust.tex",custom.note = "All models include country-of-residence and wave fixed effects. SEs are clustered by ethnic group.",caption = "Kinship Models with Trust",scalebox = .9,include.ci=F,override.pvalues = list(n1a$p.value/2,n1b$p.value/2, n1c$p.value/2, n1d$p.value/2, n1e$p.value/2),
               table = F,use.packages = F,
               custom.model.names = c("Trust", "Cultural Attitudes",  "Cultural Attitudes",  "Economic Attitudes",  "Economic Attitudes"))


texreg::screenreg(list(n1a, n1b,n1c,n1d,n1e), custom.coef.names = c("Kinship Score","Gender", "Age","Date entered in Ethnographic Atlas","Small Scale Group","Agriculture", "Malaria Index", "Trust"),custom.note = "All models include country-of-residence and wave fixed effects. SEs are clustered by ethnic group.",caption = "Kinship Models with Trust",scalebox = .9,include.ci=F,override.pvalues = list(n1a$p.value/2,n1b$p.value/2, n1c$p.value/2, n1d$p.value/2, n1e$p.value/2),
                  table = F,use.packages = F,
                  custom.model.names = c("Trust", "Cultural Attitudes",  "Cultural Attitudes",  "Economic Attitudes",  "Economic Attitudes"))

